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Bone adaptation occurs as a response to external loadings and involves bone resorption 
by osteoclasts followed by the formation of new bone by osteoblasts. It is directly trig- 
gered by the transduction phase by osteocytes embedded within the bone matrix. The 
bone remodeling process is governed by the interactions between osteoblasts and osteo- 
clasts through the expression of several autocrine and paracrine factors that control bone 
cell populations and their relative rate of differentiation and proliferation. A review of the 
literature shows that despite the progress in bone remodeling simulation using the finite 
element (FE) method, there is still a lack of predictive models that explicitly consider the 
interaction between osteoblasts and osteoclasts combined with the mechanical response 
of bone. The current study attempts to develop an FE model to describe the bone remod- 
eling process, taking into consideration the activities of osteoclasts and osteoblasts. The 
mechanical behavior of bone is described by taking into account the bone material fatigue 
damage accumulation and mineralization. A coupled strain-damage stimulus function is 
proposed, which controls the level of autocrine and paracrine factors. The cellular behav- 
ior is based on Komarova et al.'s (2003) dynamic law, which describes the autocrine and 
paracrine interactions between osteoblasts and osteoclasts and computes cell population 
dynamics and changes in bone mass at a discrete site of bone remodeling. Therefore, 
when an external mechanical stress is applied, bone formation and resorption is governed 
by cells dynamic rather than adaptive elasticity approaches. The proposed FE model has 
been implemented in the FE code Abaqus (UMAT routine). An example of human proxi- 
mal femur is investigated using the model developed. The model was able to predict final 
human proximal femur adaptation similar to the patterns observed in a human proximal 
femur. The results obtained reveal complex spatio-temporal bone adaptation. The proposed 
FEM model gives insight into how bone cells adapt their architecture to the mechanical 
and biological environment. 
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INTRODUCTION 

Bone remodeling is a dynamic process in which old bone is 
removed by osteoclasts and new bone is added by osteoblasts. 
Osteoblasts produce inorganic calcium phosphate, which is con- 
verted to hydroxyapatite, and an organic matrix consisting mainly 
of type I collagen, and then they deposit new bone on the part of 
the bone resorbed by osteoclasts. Osteoclasts dissociate calcium by 
secreting acid and degrade organic components by releasing lyso- 
somal enzymes. Interactions between osteoblasts and osteoclasts 
are critical in the regulation of bone remodeling. These coupled 
activities take place in a basic multicellular unit (BMU) (Frost, 
2001) and are modulated by mechanical and biological factors 
(Komarova et al., 2003). It is well-known that the development 
and activity of osteoclasts are under the control of the osteoblasts 
(Rodan and Martin, 1981; Hambh and Rieger, 2012). 

Disruptions in bone remodeling contribute to the pathogenesis 
of disorders such as osteoporosis and Paget's disease. Therefore, in 
order to develop more appropriate mechanobiological computer 



models to simulate the bone remodeling process, it is necessary 
to incorporate the combined effects of bone cell activities and the 
mechanical behavior of bone. 

Several categories of bone remodeling models have been pro- 
posed by different authors: models based on the global optimality 
criterion (Hollister et al., 1994; Bagge, 2000; Tovar et al, 2006; Jang 
and Kim, 2008, 2010), models based on maintaining a homeosta- 
tic state of stress/strain/strain energy (Carter et al., 1987, 1989; 
Huiskes et al., 1987; Beaupre et al., 1990; Prendergast and Tay- 
lor, 1994; Martin, 1995; Hart and Fritton, 1997; Jacobs et al, 
1997; Fernandes et al, 1999; Ruimerman et al, 2005; Hambh 
et al, 2009, 2011; Adachi et al, 2010; Hambli, 2011), models 
based on damage accumulation/repair (Prendergast and Taylor, 
1994; Martin, 1995; Ramtani and Zidi, 2001; McNamara and Pren- 
dergast, 2007), mechanistic models considering both mechanical 
and metabolic factors in the remodeling loop (Hernandez et al., 
2000, 2001; Huiskes et al, 2000; Hazelwood et al, 2001; Tay- 
lor and Lee, 2003; Taylor et al, 2004; Aznar et al., 2005), and 



www.frontiersin.org 



April2014|Volume2|Article6|1 



Hambli 



Bone remodeling based on cells activities 



a recent model considering the interstitial fluid flow (Tsubota 
et al., 2009). Some authors have developed simple 2D finite ele- 
ment (FE) remodeling simulation based on a reaction-diffusion 
system influenced by mechanical stress (Matsuura et al., 2002, 
2003; Tezuka et al, 2003, 2005). These continuum models have 
achieved some success in predicting normal bone architecture. 
They have however a major deficiency. They use mechanical stress 
or strain as a control system in which bone functional adap- 
tation is driven by the error between a mechanical set point 
and a mechanical stimulus to predict bone remodeling behavior, 
without considering the biophysical activities of osteoblasts and 
osteoclasts. 

There are currently a limited number of mathematical mod- 
els, which describe the activity of BMUs. The work by Komarova 
et al. (2003) was the first to model mathematically the non-linear 
autoregulation between osteoblasts and osteoclasts by expressing 
autocrine and paracrine factors. Rattanakul et al. (2003) proposed 
a model considering PTH as the main regulatory element in bone 
formation and resorption. In his work, Moroz et al. (2006) pro- 
posed a dynamic model that includes Michaelis-Menten type of 
feedback mechanisms. Lemaire et al. (2004) developed a more 
sophisticated model to describe the explicit molecular interac- 
tion and autoregulation between osteoblasts and osteoclasts. This 
model includes the well-known cytokine receptor activator of 
nuclear factor kB (RANK), its ligand (RANKL), and osteopro- 
tegerin pathway (OPG) (RANK/RANKL/OPG), PTH, and also 
transforming growth factor P (TGPP). Based on Lemaire et al. 
(2004), Maldonado et al. (2006) built a model, which takes the 
influence of the osteocyte under mechanical stimulation into 
account. Pivonka et al. (2008) subsequently developed a model, 
which also exhibits the RANK/RANKL/OPG pathway PTH and 
TGpp but is based on Hill functions, which are better suited 
to express the binding mechanism between ligand and receptor. 
Finally, the recent model by Ryser et al. (2009) provides enhanced 
modeling of autocrine and paracrine factors following Komarova s 
model (2003). Based on the spatio-temporal dynamic observation 
of BMU behavior, it also includes the explicit description of the 
RANK/RANKL/OPG pathway A review of the previous BMUs 
models shows that (i) the Komarova, Moroz, and Rattanakul mod- 
els are simple but need a limited number of parameters (8-2) and 
(ii) the Lemaitre, Pivonka, and Ryser models are more sophisti- 
cated but need a significantly higher number of parameters ( > 90) . 
While these cell-based models give theoretical insight into bone 
regulation mechanisms including metabolic factors, none of these 
models were implemented into an FE codes to simulate the bone 
remodeling process from a mechanobiological perspective consid- 
ering the cells activities' interactions with the mechanical reaction 
of bone. 

In current work, an extension of Komarova et al. (2003) model 
was implemented into an FE code to simulate the remodeling 
process from a mechanobiological point of view. 

The bone adaptation approach used in this study allows for 
the computation of changes in bone mass at a discrete site of 
bone remodeling at a macroscopic scale. The modeling of these 
interactions with biological factors suggested by Komarova et al. 
(2003) was completed by Bonfoh et al, 20 11) who considered the 
influence of an external loading (stimulus effects on the bone cell 



dynamics) . We combined Komarava et al.'s model and that of Bon- 
foh et al. (2011) (i) to include more general mechanical behavior 
of the bone such as fatigue damage growth and repair, mineraliza- 
tion, porosity, and bone material properties evolution and (ii) to 
include the principle of cellular accommodation, suggesting that 
the reference stimulus value for bone remodeling activation is not 
constant, but dependent on the load history. In addition, a sen- 
sitivity analysis (SA) was performed to investigate the impact of 
the model factors' sensitivities on the predicted bone density of a 
selected region of interest (ROI) (femur neck). 

The focus here was to develop and test the mechanobiological 
remodeling algorithm rather than to investigate the remodeling 
process of a real 3D proximal femur and/or develop a paramet- 
ric study of the role of the remodeling factors on bone density 
variations. The predictive potential of the current model enables 
one to investigate the effect of bone cell rate changes combined 
with mechanical external loads. Specifically, the model may offer 
a computer simulation framework to explore the development of 
new therapeutic treatments to pathological conditions and bone 
disorders such as osteoporosis. 

REMODELING MODEL DESCRIPTION 

During bone remodeling, the applied external load is transmitted 
in the form of stress/strain to the local bone site. Then the mechan- 
ical signals (stimuli) are received by osteocytes, which subsequently 
stimulate osteoclast and osteoblast populations in BMUs to change 
the bone mass. 

The corresponding mechanobiological remodeling algorithm is 
illustrated in Figure 1 . The model was implemented in the Abaqus 
code (UMAT subroutine) using a time step of 1 day. 

MECHANICAL BEHAVIOR 

To describe the continuum mechanical behavior of bone during 
the remodeling process considering the fatigue damage effects, the 
concept of continuum damage mechanics (CDM) can be used. In 
this case, the behavior law coupled to damage can be expressed by 
HambU (2010), Hambh and Thurner, 2013): 



(1) 



where aij is stress, D^^^ is the fatigue damage variable, is strain, 
and aijki is the isotropic elasticity stiffness tensor. 

For high cycle fatigue under purely elastic strain, Chaboche 
(1981) proposed a non-linear damage model given by: 



D 



fat 



1 - 



-(^) 



(2) 



where y and P are material parameters and ATf is the cycle at failure, 
which can be obtained as described by Martin et al. (1998): 

= 1.479 X 10"^^ As"^^-^ for compressive loads (3a) 
= 3.630 X 10"^^ As"^^-^ for tensile loads (3b) 

where As is the amplitude of the applied microstrain. 
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FIGURE 1 I Schematic representation of bone remodeling based on BMU every bone site. The stimulus is converted into signals, which control the 

activity coupled to mechanical stimulus: at the remodeling cycle (n), the osteoblast and osteoclast interactions. Bone formation and removal is 

applied load generates mechanical stress, strain, and fatigue damage performed by competition between osteoblast and osteoclast growth at the 

states at every FE of the mesh. A stimulus is then sensed by osteocytes at given bone site. 



Combining damage definition and elastic modulus evolution, 
the elastic modulus E at each location is calculated from the density 
p [computed using Komarovas model based on Eq. (12 and 17)], 
the bone mineralization and the damage according to Hambli et al. 
(2009): 

E=c(^l- D^^'^ p^a^ (4) 

where C is experimentally derived constants (2<p<3) and 
q = 2.74 (Hernandez et al, 2001). 

a is the ash function, which denotes the degree of bone 
mineralization expressed by Hernandez et al. (2000): 

01 (0 = amax + (OiO - a^ax) ^""^^ (3) 

0^0 5 cimax) and k denote initial mineralization, maximum degree 
of mineralization, and a parameter determining the shape of the 
temporal evolution curve. The average value for ao is about 0.65 
(Martin et al, 1998) with 0 < a < amax = 1- 

Apparent bone porosity p can be directly approximated by 
Hernandez et al. (2000, 2001): 



MECHANICAL STIMULUS 

Various expressions of the mechanical stimulus involved in bone 
remodeling have been proposed in the literature. The mechanical 



stimulus used here is expressed in terms of strain energy density. 
Therefore, the mechanical signal sensed by an osteocyte k at its 
location is given by Mullender and Huiskes (1995): 

Noc _ 

S(x, t) = ^/fc(x)^Jl^(5fc-5fc) (7) 
k=l 

where |Xk is the mechanosensitivity of the osteocyte k and Nqc 
is the number of osteocytes. Sj^ is the threshold value of the sig- 
nal considering that an equilibrium state can be obtained for near 
the reference value, 5^^. Far from this value, unbalanced activity 
of osteoblasts and osteoclasts is observed, which leads to bone 
apposition or resorption. 

fi{x) is a spatial influence function defined by Mullender and 
Huiskes (1995): 

fk(x) = exp(-dk(x)/do) (8) 

where dj^ix) is the distance between the osteocyte k and the bone 
surface location x. The parameter do is a normalization factor that 
limits the area of influence of the osteocyte. 

S]^ is the local stimulus value expressed in terms of coupled 
strain-damage energy density and is expressed by: 

Sk = l{^- D^'') cTijSij (9) 
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Osteocytes respond to low bone deformation (disuse) or 
fatigue -related microfractures by activating osteoclasts to resorb 
locally the bone matrix. This phase is followed by a period of 
bone formation by osteoblasts. Hence, mechanotransduction is 
a stimulatory signal to osteoblasts and an inhibitory signal to 
osteoclasts. 

This stimulus function coupled to damage ensures damage 
repair. When the fatigue damage reaches its critical value at fracture 
(Df^t^l)ata given bone site, the stimulus is set to zero. Therefore, 
the osteoclasts are activated to resorb the damaged bone zone. 

In addition, it is well-known that physical activity causes rela- 
tively larger changes in bone mass and strength in young people 
than in adults (Kassem et al., 1996), suggesting that the loss of bone 
in the aging may be attributed to a reduced sensitivity of bone set- 
points to the mechanical stimulus sensed by bone cells (Turner, 
1999; Frost, 2001) introduced the principle of cellular accommo- 
dation based on bone stimulus setpoint changes. He hypothesized 
that the stimulus setpoints are not constant, but are dependent on 
the load history, which can modify their values in an adaptive way 
resulting from the transient nature of many cellular biochemical 
responses to mechanical loading. Variation of the setpoints with 
time (aging) can be expressed by Schriefer et al. (2005): 



Sk = Sl + (S, - si) (1 - e-^') 



(10) 



5^ denotes the initial setpoint value and the parameter \ controls 
the velocity of the adaptation. 

BONE CELLS DYNAMIC BEHAVIOR 

Bone remodeling involves bone resorption by osteoclasts fol- 
lowed by the formation of new bone by osteoblasts. During bone 
remodeling, osteoclasts and osteoblasts interact with each other 
by expressing autocrine and paracrine factors that regulate the cell 
population. In the current work, a bone remodeling model devel- 
oped by Komarova et al. (2003) to describe the dynamics of cell 
populations at a remodeling site has been implemented into the 
FE code Abaqus. In this model, the osteoblast and osteoclast cell 
growth rates are described in the form of two differential equations 
regulated by autocrine and paracrine interactions. Autocrine sig- 
naling represents the feedback from osteoclasts and osteoblasts to 
regulate their respective formation. Paracrine signaling represents 
the factors produced by osteoclasts that regulate osteoblast forma- 
tion, and vice versa. Among the multitude of biochemical factors, 
only OPG/RANK/RANKL and TGFp/IGF pathways were modeled 
implicitly by omarova et al. (2003) in the form of non-linear 
interactions between osteoclasts and osteoblasts populations. 

The system of differential equations describing the osteo- 
clast and osteoblast rates and interactions using parameters, 
which characterize the autocrine and paracrine factors can be 
expressed by: 



dXB 

dt 

dxc _ 



gl2 gll 



dt 



B 



(11) 



where xq and xb denote, respectively, the osteoclast and osteoblast 
populations. 



ai is the osteoclast production rate, Pi is osteoclast removal rate, 
a2 is the osteoblast production rate, P2 is the osteoclast removal 
rate. 

Parameter gll describes the combined effects of all the fac- 
tors produced by osteoclasts that regulate osteoclast formation 
(osteoclast autocrine regulation). 

Parameter g22 describes the combined effects of all the fac- 
tors produced by osteoblasts to regulate osteoblast formation 
(osteoblast autocrine regulation). 

Parameter g\2 describes the combined effects of all the factors 
produced by osteoclasts that regulate osteoblast formation, such 
as TGFp (osteoclast- derived paracrine regulation). 

Parameter gll describes the combined effects of all the factors 
produced by osteoblasts that regulate osteoclast formation, such 
as OPG and RANKL (osteoblast- derived paracrine regulation). 

The model assumes that osteoclast and osteoblast apopto- 
sis is not affected by additional autocrine/paracrine regulators. 
Therefore, the decay terms of osteoclasts and osteoblasts are linear. 

The variation in bone density p at the remodeling site is 
expressed in terms of percentage of the initial mass depending 
on the number of osteoclasts and osteoblasts: 



dt 



= k2XB - hXc 



(12) 



where Ki and K2 are the normalized activities, Xc and X^ are, 
respectively, the numbers of actively resorbing osteoclasts and 
forming osteoblasts at a remodeling site defined by Komarova 
et al. (2003): 



Xc 
Xc 



xc if x^ > xc 
if Xq < Xc 



and 



Xb = x^- Xb if x^ > xb 
Xb = 0 if^B- 



(13) 



(14) 



Where xc and xb are, respectively, the number of osteoclasts 
and osteoblasts at steady state expressed by Komarova et al. (2003): 



(1-^22) 

^c = (1)^(1) 



Y 



Y 



(15) 



where 



y = gl2g21 - {I - gU) (l-g22) 



(16) 



The new density value of the bone tissue is approximated using 
the forward Euler method by: 



Pt+At = + Ap 



(17) 



During remodeling cycles, mechanical signals received by 
osteocytes stimulate pre-osteoclast and pre- osteoblast popula- 
tions, which convert the signals into autocrine and paracrine 
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factors (gij). These factors stimulate and regulate osteoclast and 
osteoblast populations, which compete in BMUs to change bone 
mass. 

In this paper, we consider the particular condition where a bone 
cell grows normally and only influences its neighbor's activity, 
but does not produce autocrine factors {gll =g22 = 0). There- 
fore, the signal received by osteoclasts and osteoblasts influences 
the autocrine and paracrine factors' productions through the 
exponents, gij by Bonfoh et al. (2011): 

^ gl2 = Ai+5ie-^i^(^'^) (18) 
g21 =A2 + 52^-^2 5(x, t) 

Al, Bi, A2, B2, Vi, and 72 are model parameters that regulate 
the production of paracrine factors S{x, t) denotes the mechanical 
stimulus function. 

The BMUs in cortical and trabecular bone have the same bio- 
logical structure but different morphological organizations. The 
remodeling or renewal of bone tissue constitutes 25% of tra- 
becular and 2-3% of cortical bone renewal each year (up to 10 
times higher for trabecular bone) (Parfitt, 1994; Rho et al, 1998). 
Therefore, to distinguish between BMU activity in cortical and 
trabecular bone, the Komarova model was applied for both bones 
but with a bone density (Eq. 12) 10 times lower for cortical bone 

(Jt^corti ^ 0.1 X kf^^i = 1, 2). 

The model parameters are summarized in Table 1. 

SIMULATION OF FEMORAL HEAD REMODELING 

To illustrate the capabilities of the mechanobiological bone adap- 
tation model developed, remodeling of a 2D proximal femur was 
performed. The 2D model is based on the geometry of a real femur, 
taken from a radiograph of a coronal section of the proximal femur 
(Jacobs et al, 1995). An FE model was constructed including the 
trabecular and cortical bones and a representation of the acetabu- 
lum allowing for the free femur head articulation (Figure 2). The 
purpose of the model is to show the adaptation of the proximal 
femur only (trabecular and cortical bone). Hence, the acetabulum 
appears as an exterior element, which participates in load bearing, 
but is not part of the adaptation process in the algorithm. 

A 2D mesh was generated using four-node plane stress ele- 
ments. Although clinical observations provided valuable informa- 
tion regarding the changes of cortical and trabecular bone types 
during growth or aging, little is known about how each of the 
cortical and trabecular bone types changes and how they interact 
with each other during the bone remodeling process (Jang and 
Kim, 2010). Thus, thresholding between cortical and cancelous 
bone is modeled by two bone regions with the same mechanical 
and BMU behavior laws but with different material properties. 

The remodeling cycle chosen here is characterized by loads rep- 
resenting a daily gait cycle (Figure 2). The daily loading history 
was simulated by three load cases consisting of joint reaction and 
abductor muscle forces similar to those proposed by Carter et al. 
(1987) for normal activity. Because it takes 3-4 months for one 
remodeling cycle to complete the sequence of bone resorption, 
formation, and mineralization (Mundy, 1999), a minimum of 6- 
8 months is required to achieve a new steady-state bone mass that is 



measurable. In the present work, the simulation was run for 5 years 
(1825 days) under normal daily loading, leading to redistribution 
of the density and remodeling parameters (baseline) condition. 

To illustrate the potential of the current integrated FE model, 
three remodeling analyses were performed corresponding to three 
cases of applied forces (Table 2): 

(a) Case NFC refers to a normal force case (standard walking gait) 
(Carter et al, 1987). 

(b) Case LFC refers to a low force case (10% of the force of case 
NFC). 

(c) Case HFC refers to a high force case (150% of the force of case 
NFC). 

A very low value of applied force was used for case LFC (the 
standard walking gait loads were reduced by 90%) to generate a 
large disuse femur region to test the potential of the algorithm 
to perform bone remodeling in sites where the level of mechani- 
cal stress is low or non-existent. Such very reduced loads concern 
patients subjected to long bed rests or can be considered as the 
effects of microgravity on bone due to long- duration space flight 
for astronauts. 

The Femur model was run in alternating load and unload 
(P=ON) increments for 1825 iterations (days) with a fixed 
number of cycles per day and orientations of forces (Table 2). 

Other inputs can influence the remodeling response (frequency, 
age, drugs, etc.). The aim of the current work was to implement 
the novel remodeling FE model and to check its validity to predict 
bone adaptation processes under different conditions. The remod- 
eling simulation can be extended by including more variables and 
inputs in order to capture complex bone behavior. 

RESULTS 

To illustrate the potential of the current mechanobiological 
remodeling model, some model factors (external load intensity 
and bone cell rates) were investigated to analyze the impact of 
these parameters on the remodeling process. The results of these 
comparative tests were analyzed in terms of the apparent density 
distribution of the proximal femur. 

The remodeling algorithm was implemented in the Abaqus 
FE code (UMAT routine) to solve the bone remodeling process, 
incorporating the mechanical and BMU behavior. The iterative 
process started from constant trabecular and cortical bone densi- 
ties and ended with variable density at the end of the remodeling 
process. 

An example of the bone adaptation sequences of a femur is 
given in Figure 3. Predicted results indicate that the bone starts 
the adaptation its density after about 3 months and undergoes 
continuous adaptation to converge to a steady state after about 4- 
5 years duration. Indeed, it can be clearly observed that after about 
4 years, no additional adaptation can be observed. 

MODEL VALIDATION 

Due to lack of specific experimental data related to human 
femur remodeling process, the complete validation of the current 
mechanobiological FE model is hard to achieve. Therefore, the 
obtained results (density distribution) were compared (i) to those 
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Table 1 | Material properties for bone used for the remodeling simulation from Mullender and Huiskes (1995), Hernandez et al. (2001), Komarova 
et al. (2003), Hambli et al. (2009). 



Parameters 




Notation 


Trabecular bone 


Cortical bone 






Initial elastic nnodulus 




Eo (MPa) 


2000 


17000 


Poisson ratio 




V 


0.3 


0.3 


Initial density 




P (g/cnn^) 


0.764 


1.4 


Density coefficient 




C (g/cnn^) 


4000 


80003 


Density exponent 




P 


3 


3 


Ash exponent 




Q 


2.74 


2.74 




Fatigue paranneter 




y 


0.2 


0.2 


Fatigue exponent 






0.4 


0.4 


MINERALIZATION PARAMETERS ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 


Initial ash fraction 




ao 


0.6 


0.6 


Maxinnunn physiological value 




O^max 


0.7 


0.7 


Velocity of the nnineralization 




k (days"^ ) 


0.0003387 


0.0003387 


STIMULUS PARAMETERS ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^HH^^^^^^^^HI^^^H 


Mechanosensitivity of the osteocyte 




|jL/^ (nnnol nnnn J"^ h"^) 


0.5 


0.5 


Osteocytes density 




A/oc (nnnn"3) 


10625 


10625 


Spatial influence factor 




do (|xnn) 


0.1 


0.1 


Acconnnnodation velocity paranneter 




X (days-^) 


0.002 


0.002 


Initial setpoint value 




5° (J m-3) 


0.0025 


0.0025 




Osteoclasts 






Osteoblasts 




Notation Trabecular 


Cortical 


Notation 


Trabecular 


Cortical 


ai (osteoclasts/day) 3 


3 


a2 (osteoblasts/day) 


4 


4 


Pi (osteoclasts/day) 0.2 


0.2 


P2 (osteoblasts/day) 


0.0017 


0.0017 


ki (osteoclasts/day) 0.24 


0.024 


k2 (osteoblasts/day) 


0.02 


0.002 


/\i 1.6 


1.6 


A^ 


-1.6 


-1.6 


ei -0.49 


-0.49 


B2 


0.6 


0.6 


Yi (g/J) 16.67 


16.67 


Y2 (g/J) 


33.37 


33.37 


xc (t=0) (osteoclasts) 15 


15 


xb {t= 0) (osteoblasts) 


1 


1 



obtained by the classic phenomenological remodeling approach of 
Huiskes (Huiskes et al., 1987; Ruimerman et al., 2005) considered 
as the gold standard for the simulation of bone remodeling (Cox 
et al, 2011) and (ii) to experimental histological proximal femur 
2D sections from the literature. 

Figure 4 shows the predicted contours of bone density at dif- 
ferent sequences for three different remodeling load amplitudes. 
Bone cell parameters were assumed to be constant (control val- 
ues). As the focus here is to assess the remodeling of the proximal 
femur only, the acetabulum was removed (post-processing) from 
the figures. 

Depending on the applied load, in the bone sites where the 
stimulus was high, the cellular activity generated bone formation 
by osteoblasts and hence, an increase in density. In the areas where 
mechanical stimulus was low, the concerned areas were under 
resorption by osteoclasts or were in a steady state. The model 
clearly indicated that in the absence of mechanical stimulus, the 
bone was not completely resorbed and reached a new steady state 
after about 60% of bone loss. 




Joint forces (femoral head) 

Single-leg stance (1) 
Abduction (2) 
Abduction (3) 
Abductor muscle forces 
Single-leg stance (1) 
Abduction (2) 
Abduction (3) 



FIGURE 2 I FE of the proximal femur and boundary conditions. 
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Table 2 | Selected load conditions that simulate walking under three different conditions. 



Remodeling load case 


Cycles/day 


Femoral head (A/) 


Orientation (°) FP 


Abductor muscle forces (N) 


Orientation (°) FP 


Low force case (LFC) (LFC = 0.1 x NFC) 


6000 


(1) 232 


24 


(1)70 


28 




2000 


(2) 116 


-15 


(2) 35 


-8 




2000 


(3) 155 


56 


(3) 47 


35 


Nornnal force case (NFC) 


6000 


(1) 2317 


24 


(1) 703 


28 




2000 


(2) 1158 


-15 


(2) 351 


-8 




2000 


(3) 1548 


56 


(3) 468 


35 


High force case (HFC) (HFC = 1.5 x NFC) 


6000 


(1) 3244 


24 


(1) 984 


28 




2000 


(2) 1621 


-15 


(2) 491 


-8 




2000 


(3)2167 


56 


(3) 655 


35 



Case NFC refers to the normal force case (standard walking gait) (Carter et al., 1987). Case LFC refers to the low force case (10% of the force of case NFC) and case 
HFC refers to the high force case (150% of the force of case NFC). The model consists of cortical (420 elements) and trabecular (1350 elements) and acetabulum 
(570 elements) regions. The femur is fixed at the bottom and loaded through the acetabulum and abductor muscle forces (details in Table 2) (Carter et al., 1987). 
Variable cortical bone thickness (dark gray), cancelous bone (light gray), and acetabulum (red). 




FIGURE 3 I Predicted bone adaptation sequences in the form of apparent bone density variation in gram per cubic centimeter. 



It can be seen (Figure 4B) that globally; the density distribu- 
tions predicted by the Huiskes model are similar to those predicted 
by current model. Nevertheless, the density levels are about 20% 
lower and the bone sites affected by remodeling predicted by the 
current model are larger than those calculated by the Huiskes 
model. These differences can be explained by the facts that: (i) 
in the early stages of the remodeling process, the current model 
starts the adaptation process by bone resorption (decrease in den- 
sity). Therefore, the final results predicted a lower density than 
that obtained with the Huiskes model, which performs the adap- 
tation process by combining formation and resorption depending 
on the stimulus level on the bone site, (ii) The model parame- 
ters governing the behavior of the osteoblasts and osteoclasts have 
never been subjected to calibration by experiments. 



All main features of femur head density distribution predicted 
by the current model are more realistic compared to those pre- 
dicted by Huiskes model. During the bone remodeling process, 
current model suggests that bone formation and resorption is 
governed by cells dynamic rather than smoothed adaptive elas- 
ticity approaches. Therefore, density heterogeneity and voids 
(Figure 4A) can be observed to correspond to the experimental 
profiles (Figure 4C). 

MODEL PARAMETERS SENSITIVITY ANALYSIS 

The application of the proposed FE remodeling model requires 
about 30 mechanobiological factors (Table 1), which depend 
among others on aging, gender, pathologies, drugs intake, etc. 
Therefore, an SA was performed to investigate the impact of these 
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Normal force case (NFC) 



High force case (HFC) 



Low force case (LFC) I Normal force case (NFC) I High force case (HF( 

(^^^g^i I Density (gW) I I^^H^^ 

■ — +8;550e-02 ml W I W . 

— +8.458e-02 Inn- 1 vppr 

-+8.367e-o2 iHC- 1 year 







1. 706e+00 
9. lOOe-02 
+9. 008e-02 
8 . 917e-02 

— +8 . 825e-02 
+ 8 . 733e-02 
+8 . 642e-02 

— +8 . 550e-02 
8 . 458e-02 
8 . 367e-02 

— +8 . 275e-02 

— +8 . 183e-02 

— +8 . 092e-02 

— +8 . OOOe-02 



lnc=2.5 years 




lnc=3.5 years 





lnc=5 years 



Current model. 





Inc=5 years 



Huiskes Model. 





Babcock's jT^BL ^BBra 

Ward' s triangle ^ ^^m - ■'■ 




Osteoporotic femur from Van 
Rietbergen et al, 2003 


Helthy femur from Jacob et al., 1 997 


Helthy femur from Van Rietbergen 
et al, 2003 



Histological section. 



FIGURE 4 I Sequences of predicted density distributions (gray level) for 
three different remodeling load levels. White regions represent low bone 
density and dark regions represent high bone density. Connparison between 
current mechanobiological and mechanical models (A) prediction based on 



the current mechanobiological model, (B) prediction based on classic 
phenomenological remodeling approach of Huiskes (Huiskes et al., 1987; 
Ruimerman et al., 2005). (C) Histological proximal femur sections from the 
literature. 



factors' sensitivities on the predicted density variation in a selected 
ROI consisting on the femur neck (Figure 5). 

Due to the relative high number of the model 30 factors, a lim- 
ited preliminary one-factor SA analysis was performed in which 
only one model parameter value was varied by a given amount 
while the other parameters were kept at their reference values 
(Hambli, 2013). 

For each parameter change, a remodeling simulation was 
performed for a duration of 5 years (1825 days) under normal 
daily loading (see Simulation of Femoral Head Remodeling for 
details) and the predicted density variation in the femur neck 
was computed for each simulation. The whole remodeling input 



parameters listed in Table 1 were each varied by —50 and +50% 
with respect to their reference values. A single run of the reference 
parameter values of the model was compared to single runs of the 
model with each parameter changed individually. SA analysis con- 
sisted of 30 total runs with a total computation time about 60 h 
on a 64 GB computer. 

DISCUSSION 

In this study, we have developed a mechanobiological FE model for 
bone remodeling, which includes a number of relevant mechan- 
ical and osteoblast/osteoclast/osteocyte processes and have used 
this model to address differences in the remodeling behavior of 



Frontiers in Bioengineering and Biotechnology | Biomechanics 



April 2014 I Volume 2 | Article 6 | 8 



Hambli 



Bone remodeling based on cells activities 




FIGURE 5 I Region of interest (ROI) where the variation of bone density 
is computed for the sensitivity analysis. 



a human proximal femur in disuse or overload and for different 
bone cell rates. A coupled strain-damage stimulus function was 
implemented, which controls the level of autocrine and paracrine 
factors. The cellular behavior is based on Komarova et al.'s 
dynamic law (2003), which describes the autocrine and paracrine 
interactions between osteoblasts and osteoclasts and computes cell 
population dynamics and changes in bone mass at a discrete site of 
bone remodeling. Therefore, when an external mechanical stress 
was applied, bone formation and resorption is governed by cells 
dynamic rather then adaptive elasticity approaches. The remod- 
eling algorithm developed has been applied for different cases 
(varying remodeling loads and osteoclast/osteoblast growth). The 
model predicted realistic and plausible results concerning the bone 
density distribution. 

All main features of femur head density distribution are realistic 
including voids representative of Ward's and Babcock's triangles. 
According to Ward's classification (Whitehouse and Dyson, 1974), 
there are different features, which characterize the human prox- 
imal femur and among them the so-called Ward's and Babcock's 
triangles. Ward's triangle represents a central area where trabecu- 
lar reinforcement is absent or with lower bone mineral densities 
than other femoral parts in the medulla and in both anterior and 
posterior walls of the neck. Babcock's triangle is located between 
the principal compressive and principal tensile trabecule. In other 
regions of the head, depending on the remodeling load amplitude, 
the density decreases gradually from the highest density region 
to the Ward's triangle region located close to the femur neutral 
axis of bending (low strain level). The current mechanobiologi- 
cal algorithm model predicts different bone density distributions 
depending on the mechanical parameters (applied external loads) 
and BMU parameters (production and removal rates) in confor- 
mity with reported clinical results. Results related to the effects of 
the BMU rates imply that the remodeling algorithm is sensitive 
to variation in the production and removal rates of BMUs. It is 
well-known that the average life span of osteoblasts (^3 months) 
exceeds the life span of osteoclasts (~2 weeks) by a factor close to 
6 (Manolagas, 2000). Therefore, the model of BMUs is most sen- 
sitive to osteoclast rates, reflecting the fact that osteoclasts are very 
active in resorbing bone with a lower population number com- 
pared to that of osteoblasts. Osteoblasts are much less active and 



require more time to form bone, as shown by the RS treatment 
simulation. 

The present model incorporates relationships between 
mechanical stimuli and the osteoclasts' and osteoblasts' produc- 
tion of autocrine and paracrine factors through the exponents, 
gl2 and g2l. These relations generated changes in BMU activity, 
which makes it possible to simulate osteoporosis treatment with 
hormones, pharmaceutical agents, exercise, and combinations of 
the three over prolonged periods of time. Such a model would be 
useful for identifying optimal treatment methodologies as well as 
changes in bone strength resulting from osteoporosis treatments. 
Therefore, one of the potential applications of the current model 
is its ability to investigate the effects of variations in BMU rates 
variation as a result of a given treatment and dose. By performing 
iterative simulations on the effects of drug treatments and doses 
on bone volume, one can predict the optimal treatment strategy 
to reduce osteoporosis and fracture risk of a specific patient. 

SENSITIVITY ANALYSIS 

The results of the SA are plotted in Figure 6, which shows the 
impact that a fixed change in each model parameter has on the 
variation of the bone density in the femur neck. 

The results of the SA underscore the influence of each para- 
meter on the resulting bone density. The results indicate that 
variations of the Young modulus and density growth factors initial 
density have a significant influence on the bone density variation 
during the remodeling process. This can be related to the fact that 
these properties explicitly affect the bone stiffness and hence the 
bone tissue deformation at a given site. Therefore, the osteocytes 
detect higher deformation signals, which are transferred ultimately 
to osteoblasts, which trigger the bone formation process (Adachi 
et al, 2010; Hambh and Rieger, 2012). Note that the SA analysis 
showed that the Poisson ratio and the ash exponent parameters 
have a negligible impact. 

Concerning the effect of the fatigue damage process, the SA 
indicated that for the current simulation boundary conditions 
(walking cycles during 5 years), the damage factors have a lim- 
ited impact on the density variation. This can be explained by 
the relative short duration of the remodeling process (5 years) 
where the level of the fatigue damage accumulation remains low 
(critical value where bone damage repair is not reached to trigger 
osteoclasts to remove the damaged bone). 

Among the mineralization process factors, the SA showed that 
the bone density is sensitive to the initial ash density followed by 
the velocity of mineralization (Eq. 5). Current model suggests that 
these factors affect the elastic modulus (Eq. 4), which controls the 
deformation level of bone. 

In the proposed model, the magnitude of the signal received by 
pre-osteoblasts and pre osteoclasts depends also on the concentra- 
tion of osteocytes (Eq. 7). This concentration can vary according 
to the age, sex, type of the bone considered, etc. Present SA results 
indicated that varying the number of osteocytes with about ±50% 
generate a slight change on the bone density variation. Neverthe- 
less, the mechanosensitivity factor, the initial set point value, and 
the accommodation velocity parameter play a significant role on 
the transduction process and hence, affect significantly the bone 
remodeling results (Adachi et al., 2010). 
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GammaZ 
B2 
A2 

Bone formation rate {k2) 
Osteoblasts apoptosis rate (Beta2) 
Osteoblasts growth rate (Alpha2) 
Gannmal 
Bl 
Al 

Bone resorption rate (kl) 
Osteoclasts apoptosis rate (Betal) 
Osteoclasts growth rate (Alphal) 



I Parameter increase (+50%) 
, Parameter decrease (-50%) 



B Initial setpolnt value 
Acconnnnodation velocity paranneter 
Spatial influence factor 
Osteocytes density 
Mechanosensitivity of the osteocyte 



velocity of the nnineralization 
maxinnum physiological value 
Initial ash fraction 



Fatigue exponent 
Fatigue paranneter 



Ash exponent 
Density exponent 
Density coefficient 
Initial density 
Poisson ratio 
Initial elastic nnodulus 



-30 



-20 



-10 



10 



20 



30 



FIGURE 6 I Relative percentage change of the predicted bone density in the femur neck. Model paranneters were varied between -50 and +50%. 
(A) Osteoblasts and osteoclasts factors, (B) transduction factors, (C) mineralization factors, (D) fatigue damage factors, (E) general (mechanical) factors. 



Concerning the osteoblasts and osteoclasts activities factors, 
current SA indicated that the bone density variation is very sen- 
sitive to the variation of all osteoblasts and osteoclasts rates. 
Nevertheless, one can notice that the remodeling process is more 
sensitive to osteoclast autocrine factors than osteoblast ones. This 
can be explained by the fact that osteoclasts are very active resorb- 
ing cells with a short lifespan, which are recruited and then 
removed rapidly. In contrast, osteoblasts are much less active with 
a longer lifespan and a lower concentration changes. Consequently, 
a greater number of osteoblasts are needed in a single bone remod- 
eling site to counterbalance the resorption of bone by osteoclasts 
in the same site. The general trend is that (i) an increase in active 
osteoblast concentration implies an increase in bone formation 
and therefore prevents bone from resorption and (ii) when vary- 
ing the cells factors by ±50%, the bone density variation is not 
symmetric, i.e., the percent of density decrease is higher to density 
increase. 

In summary, the preliminary one-factor SA showed clearly that 
variation of the remodeling factors, which can be related to aging, 
gender, pathologies, drugs intake, etc, play a significant role on the 
remodeling results in terms of bone density. Therefore, the pro- 
posed model may be applied to investigate the effects of different 
diseases and therapies on femur resistance. Therefore, one of the 
potential applications of the current model is its ability to investi- 
gate the effects of variations in BMU autocrine and paracrine rates 
as a result of a given treatment and dose, calcium-PTH regulation. 



etc. By performing iterative simulations on the effects of drug 
treatments and doses on bone volume, one can predict the opti- 
mal treatment strategy to reduce osteoporosis and fracture risk of 
a specific patient. 

The current model should be interpreted in accordance with 
the limiting assumptions contained within the model. The first 
limitation in the model was that isotropic homogeneous material 
properties were assigned as an initial condition for the remodeling 
simulation. A more realistic approach would have been to establish 
initial conditions based on subject- specific data (DEXA technique, 
QCT, etc.). In spite of these limitations, as well as the idealized 
material behavior, the predicted density distribution of the femur 
still showed many architectural features that are observed clini- 
cally. Due to lack of experimental data, the complete validation of 
the current mechanobiological FE model of bone remodeling is 
hard to achieve. The second limitation concerns the asynchronous 
activities of remodeling sites, which will be subjected to expansion 
in the near future in accordance with additional development and 
with new experimental results. Also, note that the 2D plane stress 
assumption is clearly a simplification, since it cannot completely 
represent a 3D reality. A 2D geometric model does not repre- 
sent out-of-plane properties. Using 3D FE simulation, Bitsakos 
et al. (2005) showed that muscle loads can generate significant 
changes in the periprosthetic response of bone during its remod- 
eling process in the vicinity of external loads. Nevertheless, the 
hip joint forces used in the current analyses is the greatest load 
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applied on the mediolateral plane of a femur and its magnitude 
is considerably greater than other loads. Therefore, the 2D femur 
can be considered as an acceptable representation of 3D remodel- 
ing behavior. Third, the present simulations considered the fixed 
model parameters given in Table 1. These values maybe subjected 
to change due to several factors (disease, age, drugs, gender, bone 
sites, etc.). Fourth, the present SA considered one-factor analysis 
was performed in which only one model parameter value was var- 
ied by a given amount while the other parameters were kept at their 
reference values. This simple approach showed in particular that 
the bone cells rates play significant roles on the bone adaptation 
process, which may be modulated by specific bone drugs. Never- 
theless, for future general SA analysis, it is necessary to consider 
the full factorial parameters variation simultaneously for different 
femurs geometries. 
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